home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Gold Collection
/
Software Vault - The Gold Collection (American Databankers) (1993).ISO
/
cdr11
/
ddj9304.zip
/
STEREO.ZIP
/
LISTING1.C
next >
Wrap
INI File
|
1993-01-18
|
24KB
|
634 lines
[LISTING1.C]
#include <file.h>
#include <stdio.h>
#include <unixio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define BOOLEAN int
#define TRUE 1
#define FALSE 0
#define OK 0
#define NOT_OK -1
#define WHITE 1
#define LEFT_EYE 0
#define RIGHT_EYE 1
#define INTEROCULAR_DISTANCE 2.5 /* inches */
#define PIXEL_WIDTH 0.021 /* inches */
#define ACCEPT TRUE
#define REJECT FALSE
#define SCREEN_HEIGHT 512
#define NUM_DIMENSIONS 3
#define WINDOW_BORDER_WIDTH 3
#define NUM_SHADES 256
#define MAX_LINE_LENGTH 256
#define MAX_IMAGE_SIZE 200 /* 400 x 400 samples */
#define MAX_LINE_SIZE 400 /* should be set to max_image_size
if there is enough main memory to hold the whole image. */
#define ROUND(x) ((x) > 0.0 ? (int)((x) + 0.5) : (int)((x) - 0.5))
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#define LEX TRUE
#define PC FALSE
struct point_3D
{
float x, y, z;
};
typedef struct point_3D point_3D_t;
struct point_3D_ex
{
float x, y, z;
float sha, shb; /* shading values for triangle A and B */
};
typedef struct point_3D_ex point_3D_ex_t;
extern int num_samples; /* number of samples per scan line */
extern int log_input; /* 0 if data is linear, 1 if log */
extern float z_scale; /* vertical scale */
extern float scan_sz; /* 0 < X, Y < scan_sz in nm */
extern point_3D_t min, max; /* min and max of input data */
extern point_3D_ex_t image[ MAX_IMAGE_SIZE ][ MAX_IMAGE_SIZE ];
/* Clipping rectangle boundaries - accessable to all routines */
extern double y_bottom, y_top, x_right, x_left, z_front, z_back;
extern double y_bottom_d, y_top_d, x_right_d, x_left_d, z_front_d, z_back_d;
extern point_3D_t eye_pt;
extern point_3D_t light_source;
extern point_3D_t vup;
extern point_3D_t vpn;
extern double proj_plane;
/* some global variables */
extern int f_color;
/* Local variables */
static double cos_theta, sin_theta; /* to speed up rotation, dramatically! */
#if PC
/*------------------------------------------------------------------------
LEX/90 compatible line drawing command.
--------------------------------------------------------------------------*/
void dsvec( x0, y0, x1, y1, color )
int *x0, *y0, *x1, *y1, *color;
{
if ( *color > 127 )
hline( 1, *x0, *y0, *x1, *y1 );
else
hline( 0, *x0, *y0, *x1, *y1 );
}
#endif
/*------------------------------------------------------------------------
Procedure that sets the clipping rectangle limits.
The projection plane is also set to the mid-Z coordinate.
--------------------------------------------------------------------------*/
void set_clip_volume( x_min, x_max, y_min, y_max, z_min, z_max )
double x_min, x_max, y_min, y_max, z_min, z_max;
{
y_bottom = y_min;
y_top = y_max;
x_right = x_max;
x_left = x_min;
z_front = z_min;
z_back = z_max;
}
/*------------------------------------------------------------------------
Procedure that reads the STM image data.
--------------------------------------------------------------------------*/
static void read_stm_image_data( fp )
FILE *fp;
{
register i, j;
short z_data[ MAX_LINE_SIZE ]; /* 16 bit numbers */
int count, num_lines, num_lines_loaded;
double tmp, tmp1, step, hiw; /* half image width */
long offset;
tmp = z_scale / 16385.0; /* elliminates divisions */
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
step = scan_sz / num_samples;
hiw = (double)( step * num_lines / 2.0);
offset = (long)( 2048 );
if ( fseek( fp, offset, 0 ) != 0 )
printf( "Failed to skip the STM data file header!!!\n" );
num_lines_loaded = 0;
for( i = 0; i < num_lines; i++ ) /* assumes a square image */
{
/* Read the whole line even though some of it may be thrown away. */
count = fread( (char *)( z_data ), sizeof( short ), num_samples, fp );
if ( count != num_samples ) /* reached EOF */
break;
else num_lines_loaded++;
/* Assume the data is linear, for now. */
for( j = 0; j < num_lines; j++ )
image[ i ][ j ].z = (float)(z_data[ j ]) * tmp;
#if 0
else if ( log_input == 1 )
{ /* data is logarithmic (base e) */
tmp1 = 1.0 / 10.0; /* 10 angstroms in 1 nanometer. */
for( j = 0; j < num_lines; j++ )
{
/* When STM data is logarithmic, I believe that the magnitude
is logarithmic and then the sign is added. Thus, to expontiate
properly we must expontiate the absolute value of the data. */
if ( z_data[ j ] < 0 )
image[ i ][ j ].z = (float)( -exp( (double)(-z_data[j])
* tmp ) * tmp1 );
else
image[ i ][ j ].z = (float)( exp( (double)(z_data[j])
* tmp ) * tmp1 );
}
}
#endif
/* Center the image around the origin, in the X and Y dimensions.
Setup X such that dx > 0, setup Y such that dy < 0. */
for( j = 0; j < num_lines; j++ )
{
image[ i ][ j ].x = (float)( j * step - hiw );
image[ i ][ j ].y = (float)( hiw - i * step );
}
}
max.x = min.x = image[ 0 ][ 0 ].x;
max.y = min.y = image[ 0 ][ 0 ].y;
max.z = min.z = image[ 0 ][ 0 ].z;
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
{
if ( image[ i ][ j ].x > max.x ) max.x = image[ i ][ j ].x;
if ( image[ i ][ j ].x < min.x ) min.x = image[ i ][ j ].x;
if ( image[ i ][ j ].y > max.y ) max.y = image[ i ][ j ].y;
if ( image[ i ][ j ].y < min.y ) min.y = image[ i ][ j ].y;
if ( image[ i ][ j ].z > max.z ) max.z = image[ i ][ j ].z;
if ( image[ i ][ j ].z < min.z ) min.z = image[ i ][ j ].z;
}
printf( "%f < X < %f\n", min.x, max.x );
printf( "%f < Y < %f\n", min.y, max.y );
printf( "%f < Z < %f\n", min.z, max.z );
printf( "Loaded %d scan lines of image data.\n", num_lines_loaded );
}
/* Procedure to find MIN/MAX of the image. */
find_min_max( s )
point_3D_ex_t s[][ MAX_IMAGE_SIZE ];
{
register i, j, num_lines;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
max.x = min.x = s[ 0 ][ 0 ].x;
max.y = min.y = s[ 0 ][ 0 ].y;
max.z = min.z = s[ 0 ][ 0 ].z;
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
{
if ( s[ i ][ j ].x > max.x ) max.x = s[ i ][ j ].x;
if ( s[ i ][ j ].x < min.x ) min.x = s[ i ][ j ].x;
if ( s[ i ][ j ].y > max.y ) max.y = s[ i ][ j ].y;
if ( s[ i ][ j ].y < min.y ) min.y = s[ i ][ j ].y;
if ( s[ i ][ j ].z > max.z ) max.z = s[ i ][ j ].z;
if ( s[ i ][ j ].z < min.z ) min.z = s[ i ][ j ].z;
}
}
/*------------------------------------------------------------------------
Procedure that reads the STM file header (ASCII).
--------------------------------------------------------------------------*/
static void read_stm_header( fp )
FILE *fp;
{
char line[ MAX_LINE_LENGTH + 2 ];
while( fgets( line, MAX_LINE_LENGTH, fp ) != NULL )
{
/* At any time in the header, Control-Z => end of header */
if ( strchr( line, 0x1a ) != NULL ) break;
if ( strncmp( line, "data_type", strlen("data_type")) == 0 )
{
if ( strncmp( line, "data_type = 0", strlen("data_type = 0")) != 0 )
{
printf( "File data_type is not 0 => not STM data.\n" );
printf( "Can't handle non-STM data.\n" );
exit( 1 );
}
}
if ( strncmp( line, "num_samp", strlen("num_samp")) == 0 )
{
sscanf( line, "num_samp = %d", &num_samples );
if ( num_samples > MAX_IMAGE_SIZE )
{
printf( "Warning!!! The whole image can not be processed.\n" );
printf( "The image is %d x %d.\n", num_samples, num_samples );
printf( "The maximum allowable is %d x %d.\n", MAX_IMAGE_SIZE,
MAX_IMAGE_SIZE );
printf( "Only the upper left portion of the image" );
printf( " will be processed.\n" );
}
}
if ( strncmp( line, "log_input", strlen("log_input")) == 0 )
sscanf( line, "log_input = %d", &log_input );
if ( strncmp( line, "z_scale", strlen("z_scale")) == 0 )
sscanf( line, "z_scale = %f", &z_scale );
if ( strncmp( line, "scan_sz", strlen("scan_sz")) == 0 )
sscanf( line, "scan_sz = %f", &scan_sz );
}
}
/*------------------------------------------------------------------------
Procedure that reads the STM image database.
--------------------------------------------------------------------------*/
read_stm_database()
{
char fname[40];
FILE *fp;
/* Prompt the user for the database file name */
printf( "Database file name? " );
scanf( "%s", fname );
if (( fp = fopen( fname, "r" )) == NULL )
{
printf("Couldn't open file %s.\n", fname );
return( TRUE );
}
read_stm_header( fp );
printf( "num_samples = %d, log_input = %d, z_scale = %f, scan_sz = %f\n",
num_samples, log_input, z_scale, scan_sz );
fclose( fp );
if (( fp = fopen( fname, "rb" )) == NULL )
{
printf("Couldn't open file %s.\n", fname );
return( TRUE );
}
read_stm_image_data( fp );
fclose( fp );
}
/*---------------------------------------------------------------------
Procedure to draw a border of the clipping window. The border is
completely on the inside of the window.
-----------------------------------------------------------------------*/
void draw_clip_window()
{
register i;
int x_min, x_max, y_min, y_max, color;
color = 255;
x_min = ROUND( x_left );
x_max = ROUND( x_right );
y_min = ROUND( y_bottom );
y_max = ROUND( y_top );
for( i = 0; i < WINDOW_BORDER_WIDTH; i++ )
{
dsvec( &x_min, &y_min, &x_max, &y_min, &color );
dsvec( &x_max, &y_min, &x_max, &y_max, &color );
dsvec( &x_max, &y_max, &x_min, &y_max, &color );
dsvec( &x_min, &y_max, &x_min, &y_min, &color );
x_min++; y_min++; x_max--; y_max--;
}
}
/*---------------------------------------------------------------------
Procedure that multiplies two matrices together.
Assumes that memory has been proper
-----------------------------------------------------------------------*/
matrix_mult( m1, r1, c1, m2, r2, c2, result )
double m1[];
int r1, c1; /* number of rows and columns in matrix 1 */
double m2[];
int r2, c2;
double result[]; /* resultant matrix */
{
register i, j, k;
double c;
if ( c1 != r2 ) /* Can't be done! */
return( 1 );
/* The resultant matrix is (r1 x c2). So, make every element in the
resultant matrix - one element at a time. */
for( i = 0; i < r1; i++ )
for( j = 0; j < c2; j++ )
{
for( c = 0.0, k = 0; k < c1; k++ )
c += m1[ i * c1 + k ] * m2[ k * c2 + j ];
result[ i * c2 + j ] = c;
}
return( 0 );
}
/*------------------------------------------------------------------------
Procedure that determines the average value in the image data.
--------------------------------------------------------------------------*/
float average()
{
register i, j;
int num_lines;
double sum;
sum = 0.0;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
sum += image[ i ][ j ].z;
return( (float)( sum / ( (double)num_samples * (double)num_samples )));
}
/*------------------------------------------------------------------------
Procedure that displays the vertex database as lines between successive
vertices.
--------------------------------------------------------------------------*/
display_stm_database()
{
short x, y, num_lines, intensity;
float av;
short arg[ MAX_IMAGE_SIZE * 3 ];
double z_range;
int i;
z_range = max.z - min.z;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
for( y = 0; y < num_lines; y++ )
{
for( x = i = 0; x < num_lines; x++ )
{
intensity = ROUND(( image[ y ][ x ].z - min.z ) / z_range
* (double)( NUM_SHADES ));
arg[i++] = x;
arg[i++] = y;
arg[i++] = intensity;
}
dsrnw( &num_lines, arg );
}
return( OK );
}
/*-------------------------------------------------------------------
Procedure that lets the user set the light source direction vector.
---------------------------------------------------------------------*/
void set_light_source()
{
point_3D_t ls;
double d;
printf( "Present normalized light source vector is " );
printf( "%f, %f, %f\n", light_source.x, light_source.y, light_source.z );
printf("Specify light source direction vector from point to the origin.\n");
printf( "Enter x, y & z: " );
scanf( "%f%f%f", &(ls.x), &(ls.y), &(ls.z));
if ( ls.z < 0.0 )
{
printf( "Warning! Z must be > 0\n" );
printf( "Light source direction unchanged.\n" );
return;
}
/* Make it a unit vector to allow easier interpretation of dot product
results. */
d = sqrt((double)( ls.x * ls.x + ls.y * ls.y + ls.z * ls.z ));
light_source.x = (float)( ls.x / d );
light_source.y = (float)( ls.y / d );
light_source.z = (float)( ls.z / d );
printf( "Negative normalized light source vector " );
printf( "%f, %f, %f\n", light_source.x, light_source.y, light_source.z );
}
/*-------------------------------------------------------------------
Procedure that lets the user set the viewer eye point.
---------------------------------------------------------------------*/
void set_viewer_eye_pt()
{
printf( "Present viewer eye point is " );
printf( " %f, %f, %f\n", eye_pt.x, eye_pt.y, eye_pt.z );
printf( "Specify viewer eye point.\n");
printf( "The viewing direction will be from this point to origin.\n" );
printf( "Enter x, y & z: " );
scanf( "%f%f%f", &(eye_pt.x), &(eye_pt.y), &(eye_pt.z));
printf( "Viewer eye point is" );
printf( " %f, %f, %f\n", eye_pt.x, eye_pt.y, eye_pt.z );
/* Set the view plane normal vector (from the eye_pt to origin). */
/* Vectors are defined as arrow point minus tail point. */
vpn.x = -(eye_pt.x); vpn.y = -(eye_pt.y); vpn.z = -(eye_pt.z);
printf( "View plane normal vector is" );
printf( " %f, %f, %f\n", vpn.x, vpn.y, vpn.z );
}
/*-------------------------------------------------------------------
Procedure that lets the user set the view up vector to other than default.
---------------------------------------------------------------------*/
void set_view_up_vector()
{
printf( "View plane normal vector is %f, %f, %f\n", vpn.x, vpn.y, vpn.z );
printf( "Present view up vector is " );
printf( "%f %f %f\n", vup.x, vup.y, vup.z );
printf( "Specify view UP vector from the origin to a point.\n" );
printf( "Enter x, y & z: " );
scanf( "%f%f%f", &(vup.x), &(vup.y), &(vup.z));
}
/*-------------------------------------------------------------------
Procedure that lets the user set the perspective projection plane
on the negative Z-axis to other than default.
---------------------------------------------------------------------*/
void set_per_projection_plane()
{
double e_v, e_w; /* interocular distance in the view and world
coordinate systems. */
double d;
printf( "Present perspective projection plane is %lf\n", proj_plane );
printf( "For optimum viewing distance of 18 inches from the screen\n" );
e_v = INTEROCULAR_DISTANCE / PIXEL_WIDTH;
e_w = e_v / (( x_right_d - x_left_d ) / ( x_right - x_left ));
d = 18.0 * e_w / 2.5; /* from p. 82 of Dr. Hodges' thesis */
printf( "The projection plane should be set to %lf\n", -d );
printf( "Specify new perspective projection plane on neg. Z-axis.\n" );
printf( "Enter z: " );
scanf( "%lf", &(proj_plane));
}
/*-------------------------------------------------------------------
Procedure that breaks the image into triangles and computes the dot
product between the normal of each triangle and the light source vector.
This information is necessary for subsequent shading computations.
The normal to a triangle defined by P11, P21, & P22 is P11P21 x P11P22.
The stepping distance in the X and Y dimensions is known and is constant,
as the image data has not been transformed yet (assumption).
Computations are minimized by removing redundant calculations and
relying on a uniform grid (dx == dy).
---------------------------------------------------------------------*/
void compute_normals()
{
register i, j;
point_3D_t p11_p21, p11_p22, p11_p12, cross_a, cross_b;
double d, dx; /* stepping distance (dx = dy) */
double dx_sqrd; /* dx^2 */
int num_lines;
point_3D_t na, nb; /* normals to triangle A and B */
dx = scan_sz / num_samples;
num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
num_lines--;
dx_sqrd = dx * dx;
for( i = 0; i < num_lines; i++ )
for( j = 0; j < num_lines; j++ )
{
#if DEBUG
printf( "P11 %f %f %f\n", image[i][j].x, image[i][j].y,
image[i][j].z );
printf( "P21 %f %f %f\n", image[i+1][j].x, image[i+1][j].y,
image[i+1][j].z );
printf( "P12 %f %f %f\n", image[i][j+1].x, image[i][j+1].y,
image[i][j+1].z );
printf( "P22 %f %f %f\n", image[i+1][j+1].x, image[i+1][j+1].y,
image[i+1][j+1].z );
#endif
p11_p21.z = image[ i + 1 ][ j ].z - image[ i ][ j ].z;
p11_p22.z = image[ i + 1 ][ j + 1 ].z - image[ i ][ j ].z;
p11_p12.z = image[ i ][ j + 1 ].z - image[ i ][ j ].z;
#if DEBUG
printf( "dz11_21 = %f, dz11_22 = %f, dz11_12 = %f\n",
p11_p21.z, p11_p22.z, p11_p12.z );
#endif
/* It's possible to eliminate one more multiplication in the
computations below. */
na.x = dx * ( p11_p21.z - p11_p22.z );
na.y = dx * p11_p21.z;
na.z = dx_sqrd;
#if DEBUG
printf( "Na %f %f %f\n", na.x, na.y, na.z );
#endif
nb.x = (-dx) * p11_p12.z;
nb.y = dx * ( p11_p22.z - p11_p12.z );
nb.z = dx_sqrd;
#if DEBUG
printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
#endif
/* Normalize the normal vectors, since the intensity will be
proportional to the angle between light source and the normal. */
d = sqrt((double)( na.x * na.x + na.y * na.y + na.z * na.z ));
na.x = na.x / d;
na.y = na.y / d;
na.z = na.z / d;
#if DEBUG
printf( "Na %f %f %f\n", na.x, na.y, na.z );
#endif
d = sqrt((double)( nb.x * nb.x + nb.y * nb.y + nb.z * nb.z ));
nb.x = nb.x / d;
nb.y = nb.y / d;
nb.z = nb.z / d;
#if DEBUG
printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
#endif
/* Compute the dot product between the light source vector and
the normals (== to the angle between two unit vectors ).
-1 <= cos( theta ) <= 1, which will be very useful. */
image[ i ][ j ].sha = light_source.x * na.x + light_source.y * na.y +
light_source.z * na.z;
image[ i ][ j ].shb = light_source.x * nb.x + light_source.y * nb.y +
light_source.z * nb.z;
}
}
/*-------------------------------------------------------------------
Sets a square matrix to identity.
---------------------------------------------------------------------*/
void set_to_identity( m, n )
double m[];
int n; /* n x n matrix */
{
register i, j;
for( i = 0; i < n; i++ )
for( j = 0; j < n; j++ )
if ( i == j ) m[ i * n + j ] = 1.0;
else m[ i * n + j ] = 0.0;
}
/*-------------------------------------------------------------------
Copies the source matrix to the destination matrix of the same size.
---------------------------------------------------------------------*/
void copy_matrix( s, d, r, c )
double s[], d[];
int r, c;
{
register i, j;
for( i = 0; i < r; i++ )
for( j = 0; j < c; j++ )
d[ i * c + j ] = s[ i * c + j ];
}
/* Procedure that prints out a r x c matrix. */
void show_matrix( m, r, c )
double m[];
int r, c;
{
register i, j;
printf( "\n" );
for( i = 0; i < r; i++ )
{
for( j = 0; j < c; j++ )
printf( "%lf\t", m[ i * c + j ] );
printf( "\n" );
}
}
/* Procedure to clip the ready to render polygons to the viewport.
Returns ACCEPT if the tirangle is completely inside the viewport, or
REJECT if the triangle is completely outside the viewport. */
clip_to_viewport( v, n, y )
short v[];
int n, y;
{
register i;
for( i = 0; i < n; i += 2 )
{
if (( v[i] < x_left_d ) || ( v[i] > x_right_d ))
return( REJECT );
if ( y == LEFT_EYE )
{
if (( v[i+1] < y_bottom_d ) || ( v[i+1] > y_top_d ))
return( REJECT );
}
else if ( y == RIGHT_EYE )
{
if (( v[i+1] < ( y_bottom_d + SCREEN_HEIGHT )) ||
( v[i+1] > ( y_top_d + SCREEN_HEIGHT )))
return( REJECT );
}
}
return( ACCEPT );
}
/* Procedure to set viewport limits. */
set_viewport( x_min, x_max, y_min, y_max, z_min, z_max )
double x_min, x_max, y_min, y_max, z_min, z_max;
{
x_left_d = x_min;
x_right_d = x_max;
y_bottom_d = y_min;
y_top_d = y_max;
z_front_d = z_min;
z_back_d = z_max;
}
/* Procedure to open a communication channel to the LEX/90. */
open_lex()
{
short ierr, chan;
dsopn( &ierr, &chan );
if ( ierr > 0 )
{
printf( "An error condition occured on the LEX!\n" );
exit( 1 );
}
dsclr( &(-1) );
dscsl( &2, &0, &0 );
dscer();
}
/* Procedure to configure the LEX/90 as a double buffer 8 bpp system.
The LUT is initialized as well. */
configure_lex()
{
short arg[10], read_array;
/* Configure the LEX as a two buffer system - one buffer displayable
at a time (640x512), 8-plane buffers and an 8x8 LUT. */
arg[0] = 2; arg[1] = 9;
dsesc( arg, &2, read_array, &0 );
dschan( &255, &255, &255 );
dsdisp( &0, &0, &4096 ); /* enable fill */
dsllu( &0, &0, &255, &255 ); /* ramp the LUT (256 shades of gray) */
}